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AVANT-PROPOS : Le présent document constitue le chapitre 3 du traité Modélisation numérique 
discrète des matériaux granulaires, ouvrage collectif, sous la direction de Farhang Radjaï et Frédéric 
Dubois, publié dans la collection « mécanique et ingéniérie des matériaux » aux éditions Lavoisier en 2010. 
Le titre, la pagination et l'indexation des références en sont différents. Il introduit diverses définitions et 
propriétés relatives à la mécanique des assemblages granulaires de type solide, pour lesquels un réseau de 
contacts est à même de reprendre les efforts extérieurement appliqués et de maintenir léquilibre mécanique 
sous un chargement variable, et indique comment les calculs peuvent alors être menés pour déterminer 
forces et déplacements sans faire appel à l'inertie ou à des forces dépendant du temps. Ce texte contient 
quelques références aux autres chapitres du même traité mais peut se lire indépendamment. 

1 Introduction 

Dans l'immense majorité des simulations numériques de matériaux granulaires à l'échelle 
du discret, on détermine les trajectoires d'une collection de grains en résolvant les équations 
issues du principe fondamental de la dynamique, où interviennent l'inertie et les accélérations. 
C'est ainsi que fonctionnent les méthodes de dynamique moléculaire, de dynamique des contacts 
ou pilotées par événements qui sont décrites dans les chapitres suivants de ce traité. Pourtant 
dans beaucoup de situations d'intérêt pratique, l'usage est de décrire les matériaux granulaires à 
l'échelle macroscopique par la mécanique des milieux continus solides et de traiter des problèmes 
d'évolution quasi-statique, dans lesquels l'inertie n'entre pas en compte. Le système, sous charge- 
ment variable, passe alors par une succession d'états d'équilibres. Les calculs aux éléments hnis 
dans les modèles élastoplastiques de la mécanique des sols, par exemple, sont de ce type. Par 
ailleurs, les raisonnements de changement d'échelle que l'on est tenté de bâtir pour passer du 
comportement d'un assemblage de grains avec un réseau de contact donné à une loi constitutive 
de matériau continu solide se fondent aussi sur une approche quasi-statique, dans laquelle l'ob- 
jectif est de déterminer vers quel état d'équilibre voisin de l'état initial le système est conduit 
par de petits incréments de forces appliquées, sans référence au temps physique. 

Comme nous le verrons, la possibilité même d'une telle évolution quasi-statique pose ques- 
tion, et c'est pourquoi les exemples d'utilisations de méthodes quasi-statiques sont encore rares 
dans la littérature [TJ [2j [3j |4j |5] . L'approche quasi-statique se fonde sur les matrices de raideur, 
élastiques ou élastoplastiques, dont la définition et la structure sont rappelées au § (2j où on pré- 
sente également d'autres notions fondamentales de la mécanique des réseaux de contact, comme 
les matrices de rigidité (à ne pas confondre avec les matrices de raideur). On montre ensuite 
(Éj3|) comment ces différentes matrices et leurs propriétés conditionnent les configurations des 
assemblages granulaires en équilibre stable, et on rappelle des observations issues des simula- 
tions de systèmes de disques ou de sphères qui illustrent leur importance pratique et théorique. 
L'application des méthodes matricielles quasi-statiques à l'élasticité des assemblages granulaires, 
puis à leurs déformations inélastiques fait l'objet des §[J]et[5] La conclusion [6] évoque brièvement 
quelques directions nouvelles vers lesquelles on pourrait appliquer ou étendre la méthode. 
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2 Statique et cinématique des réseaux de contact 



On considère ici des grains en interaction par des contacts ponctuels ou de très faible 
étendue (on supposera par exemple que la forme des grains est régulière et strictement convexe, 
en excluant les contacts par des faces ou des arêtes de polyèdre). Les modèles mécaniques de 
matériaux granulaires utilisent en général des lois de contact, qui mettent en correspondance 
certains mouvements relatifs des grains avec les efforts de contact. Ces correspondances peuvent 
être linéarisées pour de petits incréments, de même que la cinématique, et on fait ainsi apparaître 
certaines matrices, dont le rôle est central dans l'approche quasi-statique. Pour simplifier, on ne 
donnera leur écriture complète que dans le cas de grains circulaires à deux dimensions (2D) ou 
sphériques à trois dimensions (3D). La définition de la matrice de rigidité est ici conforme à celle 
de la théorie de la rigidité (pour les treillis ou les systèmes de tenségrité [6]), que nous étendons 
aux assemblages granulaires (ce n'est pas une matrice de raideur). 

2.1 Hypothèse des petites perturbations (HPP) 




I I 

FIGURE 1 - Grandeurs associées à l'interaction entre deux grains i et j. On définit le point de 
contact, ou, pour chacun des deux grains, le point de sa surface le plus proche de son vis-à-vis. Ces 
points définissent les extrémités des vecteurs-branches Ry et Rjj, dont les origines respectives 
sont les centres (arbitrairement choisis) de i et de j. Le vecteur unitaire tijj pointe de i vers j 
et est normal aux deux surfaces lorsqu'elles sont en contact, hij désigne l'interstice ou distance 
minimale entre les deux surfaces. 

Les approches quasi-statiques sont adaptées aux calculs incrémentaux, dans lesquels on 
cherche à relier de petits incréments de forces appliquées aux petits déplacements des grains, 
dont la cinématique est celle d'une collection d'objets solides indéformables. Sauf indication 
contraire on aura recours dans la suite à V hypothèse des petites perturbations (HPP), c'est-à-dire 
que l'on fera l'approximation qui consiste à négliger l'effet des déplacements sur la géométrie de 
l'assemblage granulaire. On traite alors les déplacements comme des vitesses, et les grandeurs 
géométriques, dont la définition est rappelée sur la figure Q] sont gardées constantes. Ainsi, les 
vecteurs Ry, Rjj, restent fixes, tandis que hij sera une fonction affine des déplacements 
et des (petites) rotations des grains. Les effets de l'HPP seront discutés a posteriori. Cette 
approximation n'est pas une vraie limitation à l'approche quasi-statique, qui peut s'en dispenser. 
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Mais on verra que les erreurs qu'elle introduit sont faibles, et nous l'adoptons aussi parce que 
l'exposé des méthodes de calcul en est simplifié. 



2.2 Déplacements, déplacements relatifs et matrice de rigidité 

Nous considérons, en dimension D égale à deux ou trois, une collection de N grains. 
Chacun d'entre eux possède n; = D{D + l)/2 degrés de liberté (translation et rotation). Nous 
nous plaçons au voisinage d'une configuration de référence, à partir de laquelle les déplacements 
sont traités selon l'HPP. Les conditions aux limites peuvent faire intervenir des objets particuliers, 
comme des parois, dont certains degrés de liberté sont figés ou bien astreints à un mouvement 
extérieurement imposé. Pour simplifier, on admettra que les actions extérieures sur le système 
consistent soit à interdire certains mouvements (cas d'une paroi fixe par exemple), soit à imposer 
des forces (ou des moments) sur certains objets. On notera n g le nombre de degrés de liberté 
associés aux conditions aux limites imposées au système. Le nombre total de degrés de liberté 
est Ni = ni x N + n g , soit 6iV + n g en 3D et 3N + n g en 2D. Un exemple simple est illustré par 
la figure |2] : on considère un assemblage 2D de N grains, enfermé dans une cellule rectangulaire 
constituée de 4 parois, dont 2 (marquées 3 et 4 sur la figure) sont fixes, tandis que celles qui leur 
sont opposées (respectivement : 1 et 2) possèdent un seul degré de liberté de translation dans 
la direction qui leur est orthogonale. Cette configuration permet ainsi la compression biaxiale. 
Un cas intéressant (objet du chapitre 6 du présent ouvrage) est celui des conditions aux limites 




FIGURE 2 - Un choix de conditions aux limites en 2D adaptées à la compression biaxiale, avec 



2 degrés de liberté associés aux parois. 



périodiques. On peut alors avoir des déformations globales de la cellule de simulation, qui se 
superposent aux mouvements des grains à l'intérieur d'une cellule périodique de forme et de 
taille fixes. Ainsi, au lieu de simuler la compression biaxiale d'un échantillon 2D en le confinant 
entre des parois comme sur la figure [2j on peut se donner une cellule rectangulaire périodique, et 
écrire le déplacement Uj de chaque grain i, dont la position du centre, par rapport à une origine 
quelconque (mais qu'il est commode de placer au centre de symétrie de la cellule), est repérée 
par le vecteur rj, sous la forme : 

Uj = ûj-e-rj. (1) 
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Dans ([T]), le tenseur de déformation e (ici défini selon la convention de la mécanique des sols, 
c'est-à-dire que les raccourcissements sont positifs) a la forme diagonale 



e = 









SL 2 
L 2 



(2) 



avec de petites variations ALi, AL2 des dimensions de la boîte, traitées comme infinitésimales 
dans le cadre HPP, tandis que ûj désigne un déplacement supplémentaire qui satisfait aux condi- 
tions de périodicité. AL\ et AL2 sont les n g = 2 degrés de liberté associés aux conditions aux 
limites dans ce cas. 

La mécanique de l'assemblage granulaire est sensible aux déplacements relatifs dans les 
contacts. Comme on considère de petits mouvements dans le voisinage immédiat d'une configura- 
tion donnée, on peut se donner une liste a priori de paires de grain en contact, ou éventuellement 
susceptibles d'entrer en contact. Pour chaque paire i, j, on choisit arbitrairement le grain origine 
i et le grain extrémité j (comme quand on oriente les connexions sur un graphe), et le déplace- 
ment relatif Uij se définit comme la différence entre les déplacements du point de contact selon 
le mouvement de i et selon le mouvement de j. Si on désigne les (petites) rotations par 9i, 9j, 
les déplacements des centres étant u,, Uj, on a alors 



formule générale faisant apparaître les vecteurs-branches définis sur la figure [TJ ce qui donne 
dans le cas des grains sphériques ou circulaires de rayons Ri et Rj, en choisissant le centre 
conventionnel au centre géométrique du grain, 



Dans le cas de conditions aux limites périodiques on fera apparaître les déplacements ûj au second 
membre de ([3]) ou (j4|), ainsi qu'un terme supplémentaire impliquant les déformations globales au 
niveau de la cellule de simulation. Dans le petit exemple précédent de la compression biaxiale 
2D, on aura ainsi, pour des grains circulaires (le vecteur rotation devenant scalaire) 



où r ij désigne la plus petite image de rj — par l'une des translations qui à la cellule associent la 
famille de ses copies périodiques, la matrice e étant définie ainsi qu'en ([2j), et le vecteur unitaire 
tangentiel tij complétant pour former une base directe. 

Il est commode de définir un unique vecteur déplacement U avec autant de coordonnées 
que de degrés de liberté, iVj, en agrégeant les coordonnées des déplacements et des vecteurs 
rotations de tous les grains, de 1 à N, puis les n g degrés de liberté associés aux conditions aux 
limites (déplacements de parois, paramètres de déformation d'une cellule périodique). Avec N c 
contacts, on définit de même un vecteur des déplacements relatifs U, avec DN C coordonnées 
en dimension D. Il est d'usage d'utiliser un système de coordonnées dans lequel on isole, pour 
chaque contact, la composante normale du déplacement relatif, soit ïXij-Uij. Les coordonnées dans 
l'espace des IÀ sont donc, étant donnée a priori une liste ordonnée des contacts, le déplacement 
relatif normal, puis la ou les {D — 1) coordonnée(s) du déplacement relatif tangentiel dans le 
premier contact, puis le déplacement relatif normal dans le deuxième, etc., pour terminer par les 
incréments des n g degrés de liberté de la cellule de simulation. 

On voit alors que les relations (J3j) définissent une application linéaire U 1— > U, 



— Uj + 9. t x Rjj — Uj — 9j x Rjj, 



(3) 




(4) 




(5) 



U = G -U 



(6) 
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La matrice G est la matrice de rigidité associée à l'ensemble de grains et de contacts - la struc- 
ture granulaire (au sens de la mécanique des structures). G est une matrice à DN C lignes et Ni 
colonnes. Le noyau de G est l'espace des vecteurs déplacements U qui ne conduisent à aucun 
déplacement relatif. De tels déplacements sont aussi appelés des mouvements de mécanismes. La 
dimension de l'espace qu'ils forment, notée k dans la suite, est par définition le degré d'hypo- 
staticité de la structure. L'image de G est le sous-espace de lR DNc constitué des déplacements 
relatifs pour lesquels il est effectivement possible de trouver des valeurs de déplacements (et de 
rotations, etc.) qui leur correspondent. C'est l'espace des déplacements relatifs compatibles, de 
dimension Ni — k. 

2.3 Efforts extérieurs, forces de contact et matrice de rigidité. 

Sur chacun des grains, dont la cinématique est celle d'un objet rigide, on peut appliquer 
une force et un moment. De plus, à chacun des n g degrés de liberté associés aux parois ou à la 
cellule de simulation, soit X, on peut aussi faire correspondre une force extérieure généralisée 
J-, telle que son travail dans un « déplacement » AX soit JAX. Dans le cas du mouvement 
d'une paroi, comme sur la figure [2] il s'agit d'une force au sens le plus ordinaire. En général, on 
peut avoir affaire à des forces en un sens généralisé. Ainsi lorsque le degré de liberté de la cellule 
a le sens d'une coordonnée du tenseur des déformations, e a ^, alors la force généralisée est le 
produit du volume par la coordonnée correspondante du tenseur des contraintes, a a p. Ces efforts 
extérieurs définissent donc un vecteur F ext avec Ni coordonnées, le travail s'écrivant F ext • U. Ils 
doivent être équilibrés par les forces aux contacts qui définissent un vecteur f dans un espace de 
dimension DN C . Par définition, la force de contact F™ est la force transmise par le grain origine, 
i, au grain extrémité, j, au point de contact. L'équilibre des forces généralisées correspondant à 
tous les degrés de liberté prend la forme d'un ensemble de relations linéaires : 

F ext = g f ( ? ) 

définissant une matrice HàJVj lignes et DN C colonnes. Le théorème des travaux virtuels, dont 
la vérification directe est facile dans les cas qui nous intéressent ici, énonce que H n'est autre 
que la matrice transposée de G : 

_ H= T £. (8) 

La même matrice de rigidité apparaît donc dans les relations statiques et cinématiques. Par 
exemple, dans le cas de la compression biaxiale d'un système de disques avec des conditions aux 
limites périodiques, évoqué plus haut, les forces généralisées correspondant aux paramètres de 
déformation e a , a = 1, 2 sont Ao~ aa (a = 1, 2), A désignant l'aire du système dans la configuration 
de référence. On sait que l'on a, la somme étant étendue à tous les contacts (identifiés par la 
paire ordonnée de disques concernés) : 

Aa a0l = Y j rt j fi i r (9) 

i<j 

Cette relation est l'équation du système linéaire ([7J relative relative à force conjuguée de la défor- 
mation e a , et les coefficients de la ligne correspondante de la matrice H sont les rfj. D'après ([5j) 
ce sont aussi les coefficients de la colonne de G correspondante. 

Le noyau de H est constitué des forces de contact f autoéquilibrées, c'est un sous-espace 
de ~S{P Nc dont la dimension h est par définition le degré d'hyperstaticité de la structure (ou degré 
d'indétermination des forces). Les forces intérieures équilibrant le chargement F ext appliqué, si 
elles existent, forment un espace affine de dimension h, car à toute solution particulière de ([7]) 
on peut ajouter une solution du système homogène associé, c'est-à-dire un vecteur de forces 
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de contact autoéquilibrées. Quant à l'image de H, c'est le sous-espace de M Nl constitué des 
chargements supportables, qu'il est possible d'équilibrer avec des forces intérieures. 

L'image d'un opérateur étant égale à l'orthogonal du noyau de son transposé, on dis- 
pose d'une caractérisation des déplacements relatifs compatibles comme orthogonaux à tous les 
systèmes de forces de contacts autoéquilibrées, ainsi que d'une caractérisation des chargements 
supportables comme orthogonaux aux mouvements de mécanisme (c'est-à-dire qu'ils ne tra- 
vaillent pas dans ces mouvements). De plus on obtient une relation entre degrés d'hypostaticité 
et d'hyperstaticité : 

Ni + h = DN C + k. (10) 



2.4 Lois de contact et matrices de raideur 

Les lois de contact, dont nous supposons, conformément aux modèles les plus courants, 
qu'elles combinent élasticité et frottement de Coulomb, peuvent être mises sous une forme incré- 
mentale. Elles font alors apparaître des paramètres de raideur. A chaque contact i,j est associée 
une matrice de raideurs locales KL. ., qui relie les incréments des coordonnées de la force de contact 
aux variations du vecteur déplacement relatif : 



Af, 



Les lois de contact habituelles conduisent à décomposer fjj en ses composantes normale et tan- 
gentielle, selon 



ij ïi-ij ~\- T ij 



avec 



T« 1 n 



(12) 



Dans (|lip . on choisit par conséquent pour les deux membres des axes de coordonnées selon la 
direction normale pour le premier et selon celle de la force de contact tangentielle actuelle (avant 
incrémentation) Tjj pour le deuxième, c'est-à-dire les 3 vecteurs de base suivants : 



rijj x tij. 



(13) 



- y 



Dans les modèles simples où l'élasticité du contact est prise linéaire et unilatérale, chaque contact 
est doté de raideurs K 7 ^, Kj, indépendantes des forces qu'ils transmettent et on a pour /Ç„ la 
forme diagonale simple : 



K, 1 



=i.i 



K %3 
N 















(14) 



Ï4]) ne convient que si l'inégalité de Coulomb est satisfaite 



l'usage de la forme élastique ([TT 
sous forme stricte : ||Tjj|| < fiNij. S'il y a égalité dans la condition de Coulomb, alors les raideurs 
locales dépendront de la direction de l'incrément AUij, car il faut alors discuter selon le statut 
du contact [7J. On a donc : (tjj est défini en H3|) 








K% 



si Kj,AUij ■ — [iK^r^ij ■ AUij > 



dans le cas contraire 



(15) 



La prise en compte de l'élasticité de Hertz-Mindlin au contact [8l |9] introduit plusieurs modi- 
fications, car les raideurs Kn et Kt dépendent de la déflexion au contact hij (ou de la force 
normale) [8j[9], ainsi qu'en général de l'histoire du mouvement relatif tangentiel au contact. En 
fait la matrice des raideurs locales prend différentes formes suivant la direction de AUij même 



si le frottement n'est mobilisé nulle part dans la région du contact [9] (c'est-à-dire même si le 
coefficient de frottement \x est infini). Cette dépendance directionnelle des raideurs locales n'a 
toutefois qu'une influence très faible dans le calcul des propriétés élastiques macroscopiques, 
comme il est montré dans la référence [7], où on trouvera davantage de détails sur la forme 
des matrices de raideur avec un modèle de Hertz-Mindlin. Dans la suite on se limitera dans les 
exemples à la discussion des propriétés élastiques dans le cas de grains sphériques avec élasticité 
de Hertz-Mindlin (§ [4j ou des propriétés élastoplastiques (§ [5]) avec des disques (2D) et une 
élasticité linéaire unilatérale dans les contacts. 

Quelle que soit la forme des matrices carrées /C. (de taille D x D) associées à chacun des 
contacts, on peut les rassembler dans une grande matrice carrée de dimension DN C x DN C , la 
matrice des raideurs de contact /C, qui est diagonale par blocs, car elle ne couple pas les contacts 
différents, et qui relie des vecteurs de TR DNc , incréments de déplacements relatifs et de forces de 
contact : 

Af = g-AÛ. (16) 

Rappelons qu'en général, la relation (|16p n'est linéaire qu'en apparence puisque la matrice prend 
des formes différentes selon la direction de AU. Rassemblant les équations ([6]), (116|) . puis ([7} et 
(|8j), on obtient la relation entre incréments de déplacements et de chargement extérieur, qui fait 
apparaître la matrice de raideur IC 1 ' : 

AF ext = g(i) . AU) ayec = T g . £ . £ (17) 

À la différence de la matrice de rigidité, la matrice de raideur K^ 1 ) est carrée, avec autant de lignes 
et de colonnes que de degrés de liberté. Dans le cas où chaque bloc diagonal de JÇ est de la forme 
élastique (|14p . avec des raideurs toutes strictement positives, alors (j 1 7[) montre immédiatement 
que K^ 1 ) est symétrique et positive, et que son noyau coïncide avec celui de la matrice de rigidité 

2.5 Contributions géométriques à la matrice de raideur 

En ignorant les variations des directions normales et des vecteurs-branches dans les dé- 
placements, on néglige, dans le cadre de l'HPP, une autre contribution à la matrice de raideur 
- notons-là K^ 2 ^ - dont l'origine est la suivante. Nous étudions l'effet de petits déplacements 
au voisinage d'une configuration donnée, dans laquelle l'assemblage est déjà soumis à des forces 
extérieures. Chaque contact i-j transmet alors une force fy, qui varie lorsque les grains se dé- 
placent, d'abord en raison du mouvement relatif au point de contact, d'où l'incrément de force 
calculé avec le comportement du contact, qui s'exprime avec la matrice de raideur K^ 1 ^ ; mais 
aussi simplement parce que la force de contact doit suivre les grains dans leur mouvement, ce 
qui donne une seconde contribution à l'incrément Afjj et fait apparaître le second terme K 1 - 2 ) 
de la matrice de raideur. Si, par exemple, i et j se déplacent ensemble comme un seul corps 
rigide, il n'y a aucun mouvement relatif au contact. La force fjj doit alors suivre ce déplacement 
matériel. Si le mouvement relatif de i et j est un roulement sans glissement, c'est-à-dire une 
rotation par rapport à un axe orthogonal au vecteur normal n^-, ou un pivotement, c'est-à-dire 
une rotation autour de ny, alors aucune variation de fjj ne provient non plus du déplacement 
relatif, qui reste inchangé. La variation de la force fy dans de tels mouvements n'est curieuse- 
ment presque jamais discutée dans la littérature où sont présentées les lois de contact et la mise 
en œuvre des méthodes aux éléments discrets. Il est en général supposé que la force de contact 
suit le mouvement de rotation du vecteur unitaire normal n^- dans le roulement. Cette règle est 
automatique pour la composante normale Nijïiij, dont l'intensité Nij est généralement fonction 
de la déflexion normale hy. On doit tenir compte du terme A^jArijj dû à la rotation du vecteur 
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unitaire normal dans le calcul de l'incrément de la force de contact. Pour deux sphères centrées 
en r-j et Tj, avec r. 



et Tïj 



on a : 



Arijj = — (l - n*j (g) rijj) • (Au,- - Au,) . 



(18) 



Pour la composante tangentielle T,,-, le traitement habituel des calculs 2D, dans lequel on 
construit un vecteur unitaire tangentiel, est similaire à celui de la composante normale. Dans les 
calculs 3D, il faut explicitement faire tourner la force tangentielle avec ny, et aussi prendre en 
compte le pivotement. Pour ce faire, un choix naturel est d'imprimer à Ty la rotation moyenne 
des objets i et j autour de ny. Au total |10[ [7] cela conduit à incrémenter la force tangentielle 

(2) 

de AT\y défini pour deux grains sphériques par la règle suivante : 



AT (2) 



- y 



(Au; 



A9i + Adi • n 



(rijj x Ty) 



(19) 



Cette contribution AT^ est à ajouter au terme AT^ , partie tangentielle de Afjj calculée 
selon (|16p . En rassemblant les contributions des composantes normales et tangentielles on trouve 
pour les incréments de force de contact dus à la précontrainte 

Af {2) = g ■ AU, (20) 

où £ est une matrice à DN C lignes et Ni colonnes, dont la ligne de blocs D x ni relative au 
contact i-j ne contient que les deux éléments non nuls et Ç . ., qui s'écrivent, si on prend la 
base (jT3j) pour écrire les coordonnées de Afy et aussi de AU, et en notant Ty = ||Ty||, 



£... 
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(21) 



Pour obtenir la partie géométrique K^ 2 -* de la matrice de raideur K, il faut revenir à un système 
de coordonnées dans une base fixe, le même pour tous les blocs de £, c'est-à-dire que chacun des 
deux blocs 3x3 des matrices £ ... et £ . . . écrites en ( 1211) , le bloc correspondant aux déplacements 
comme le bloc correspondant aux rotations, est à multiplier à droite par la matrice 3x3 dont 
les vecteurs-colonnes sont ceux de (fTH]) . Une fois £ ainsi transformée on peut alors écrire 

K = + K (2 \ avec donné par (HZ]) et ^ = T £ ■ g. (22) 

C'est en toute rigueur avec la matrice de raideur complète K qu'il faut écrire la relation entre 
incrément de chargement et variation de déplacement plutôt qu'avec K^ 1 ^ seulement 



AF' 



ext 



K AU. 



(23) 



La distance r« présente dans les formules (12111 est confondue en général avec la somme des rayons 



Ri + Rj, en négligeant hij. Comme elle provient de la variation du vecteur normal selon (j!8p . 
il s'agit des rayons de courbure des surfaces en contact et non des vecteurs-branches qui, eux, 
interviennent dans la matrice de rigidité selon fl5J) . Des formules générales pour les coefficients des 
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matrices K^ 2 ^ dans le cas de grains de forme quelconque (mais régulière et strictement convexe) 
sont fournies dans les références [TT] et |12j . On notera que la matrice K*- 2 ^ n'est pas symétrique. 
Elle le devient toutefois dans le cas particulier de grains sphériques ou circulaires sans frottement, 
pour lesquels on peut ignorer toutes les rotations (qui sont autant de mouvements de mécanisme), 
car on a alors pour les blocs D x D couplant forces et déplacements : 

T%3 (24) 
K^ 2 ** = — K( 2 ^ (somme sur tous les j en contact avec i) 

j 

D'après (|2ip et (|22p . et compte tenu des coefficients de G, les coefficients de la matrice K^ 2 ) 
couplant les forces aux déplacements sont d'ordre F/R, où F est une force de contact typique 
et R un rayon de courbure ou une longueur de vecteur-branche, alors que les coefficients de 
K^ 1 ) correspondants sont les raideurs Kjy, Kt, multipliées par des coefficients d'ordre 1 (comme 
les coordonnées des n*j). Les coefficients qui couplent les forces aux rotations font apparaître un 
vecteur-branche supplémentaire et sont de l'ordre de F (et ceux de KW d'ordre K N ), et ceux des 
lignes relatives aux moments comportent également un facteur R supplémentaire. Chacun des 
coefficients de K^ 2 ^ se compare donc à son analogue dans K^ 1 ^ comme F/R aux raideurs Kn ou 
Kt- Or F est d'ordre K^h pour une déflexion de contact typique h, et on a h -C R\H d'où <C 
pour chacun des coefficients de la matrice couplant deux objets en contact I12 | 17]. Il est 
donc légitime de négliger la contribution géométrique à la matrice de raideur, et d'approximer fj23[) 
par (|17p . sauf pour les vecteurs U tels que • U = 0. 



2.6 Stabilité 

Un critère classique de stabilité d'un état d'équilibre |TT] est que la forme quadratique 
AV(AU), définie par 

A 2 W(AV) = AU • g • AU (25) 

soit positive. Comme dans (125 1> . K • AU est un incrément de chargement AF ext , on peut écrire 
la condition comme AF ext • AU > pour tout AU, c'est-à-dire que le « travail du second 
ordre » doit être positif. Pour comprendre l'origine de ce critère, notons que la relation (123p . que 
nous avons écrite en supposant l'équilibre de forces, exprime en fait que l'incrément des forces 
intérieures AF mt est égal à — K • AU. Lorsque, partant d'un état d'équilibre, une perturbation 
extérieure imprime aux grains des vitesses V (rassemblées dans un vecteur à Ni coordonnées, 
comme les petits déplacements), à temps court les grains se sont déplacés de U = Vt + 0(t 3 ) 
(puisqu'à U = on a une configuration d'équilibre par hypothèse) et subissent des forces 
intérieures AF mt = — K • Vt. La puissance AF mt • V en sera strictement négative si A 2 W est 
définie positive, d'où une décroissance de l'énergie cinétique produite par la perturbation, qui 
aura tendance à augmenter, en revanche, pour un vecteur V, s'il existe, tel que A 2 W(V) < 0. 
Une telle situation est similaire à celle d'un ressort dont la raideur serait négative, et dont la 
réponse tendrait à augmenter l'élongation plutôt que de s'y opposer par une force de rappel. 

Si K est une matrice symétrique alors on peut, pour de petits déplacements U au départ 
d'une configuration équilibrée, décrire les incréments de forces intérieures — K-U comme dérivant 
de l'énergie potentielle quadratique A 2 W (la symétrie de K assure l'égalité des dérivées secondes 
croisées). L'état d'équilibre est stable si l'énergie potentielle est minimale, ce qui requiert la 
positivité de A 2 W. Pour les matériaux granulaires la symétrie de K, et donc la définition d'une 

1. Au chapitre 9 de ce même traité on définit un paramètre sans dimension k caractérisant le niveau 
de raideur de l'assemblage granulaire, tel que h/R soit d'ordre k -1 , à partir de la pression de confinement, du 
diamètre des grains et des raideurs de contact. On aura donc K ( - 2) /K (1) = 0(k~ 1 ) en général. 
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énergie potentielle au voisinage d'un état d'équilibre sous un chargement donné, signifie qu'il 
existe, au moins pour de faibles déplacements, un régime de comportement élastique. 

3 Illustrations et discussion 

Après le rappel de diverses notions liés aux réseaux de contact nous évoquons dans cette 
section divers résultats et observations numériques qui montrent leur importance comme outils, 
pratiques (critères d'équilibre) ou théoriques (formulation de problèmes de calcul à la rupture, 
évaluation des influences de la géométrie du réseau, de la forme des grains) d'analyse et de 
compréhension de la mécanique des structures de contacts intergranulaires. 

3.1 Définitions, rôles des matrices de rigidité et de raideur 

L'approche quasi-statique étant relativement peu répandue, les définitions et la termino- 
logie ne sont pas fixées de façon unique. Ainsi les références [5j [Ï3] utilisent une matrice (dite 
de contact ou de configuration) notée c qui est définie comme — T G ici, et c'est à — T G que les 
auteurs de |T3] donnent le nom de matrice de rigidité. Quant à la matrice de raideur K elle est 
parfois appelée « matrice dynamique » |15[ 116]. comme en physique du solide. 

La matrice de rigidité G est une donnée fondamentale de la structure granulaire, elle 
apparaît naturellement dans la description de la méthode de dynamique des contacts, où elle se 
combine à une matrice d'inertie, alors qu'elle se combine à l'élasticité du contact en dynamique 
moléculaire (chapitre 2 de cet ouvrage) comme dans l'approche quasi-statique décrite ici. Elle ne 
dépend, par (|3|), que des positions des centres des grains et des points de contact. Sa transposée 
H = T G exprime l'équilibre des forces de contact par J7J). De telles forces de contact f sont dites 
statiquement admissibles. On dit que f est plastiquement admissible si l'inégalité de Coulomb 
est satisfaite dans chaque contact. On peut alors, sans se préoccuper de la forme précise de la 
loi de contact et de petites déformations éventuelles, s'intéresser à l'ensemble S des forces de 
contact à la fois statiquement et plastiquement admissibles. S est l'intersection de l'espace affine 
de dimension h (le degré d'hyperstaticité) des forces de contact statiquement admissibles et du 
cône des forces de contact plastiquement admissibles, c'est un ensemble convexe. L'approche du 
calcul à la rupture consiste à déclarer chargement supportable tout F ext pour lequel S est non 
vide. Nous verrons que cette condition, certes nécessaire, ne garantit pas la stabilité d'un réseau 
de contact sous le chargement considéré, et qu'il peut y avoir rupture avec 5^0. 

Nous admettons dans ce chapitre que les efforts de contact sont des forces ponctuelles. Il 
arrive toutefois que, pour modéliser les contacts entre grains par plusieurs aspérités dans le cas de 
surfaces rugueuses, on introduise une résistance au roulement dans les contacts |17| 118] . La loi de 
contact correspondante, reliant un moment à une rotation relative peut être prise analogue à la 
loi tangentielle |19j . avec un coefficient de frottement de roulement [Xr (une longueur) qui limite 
la valeur du moment Y à urN , et une raideur en rotation. On doit alors étendre la définition 
de f, dont la dimension passe de DN C à D{D + \)N c /2 pour y inclure des moments au contact 
(avec deux moments de roulement et un moment de pivotement en dimension 3), tandis que IÀ 
contiendra des rotations relatives. La référence |20| contient une brève discussion des matrices 
de rigidité et de raideur dans un assemblage granulaire 2D avec résistance à la rotation. 

L'écriture de la matrice de raideur sous la forme K = T G • JÇ • G + T G • Ç permet de 
dégager le rôle et l'influence des différentes données géométriques et mécaniques : la structure du 
réseau de contacts détermine la matrice de rigidité G, les lois de contact, une fois linéarisés pour 
de petits incréments, fournissent la matrice (diagonale par blocs associés à chacun des contacts) 
des raideurs de contact /C, et la forme des grains, plus précisément leur courbure aux points de 
contact, n'apparaît que dans C. Cette contribution géométrique à la matrice de raideur intervient 
dans les questions de stabilité. Enfin, le terme de pivotement au second membre de (jl9[) affecte 
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aussi la matrice C, mais son origine n'est pas liée à la courbure des surfaces. Il est nécessaire 
(bien qu'en général oublié!) pour assurer l'objectivité [llj du modèle mécanique du contact : si 
les deux grains en contact sont animés d'une même rotation de corps rigide autour du vecteur 
normal, la force tangentielle doit subir cette rotation elle aussi. 

3.2 Stabilité et équilibre 

On est très souvent confronté, dans les calculs par éléments discrets à la question de la 
tolérance avec laquelle des équations d'équilibre sont satisfaites. Lors de la simulation d'un essai 
biaxial, comme schématisé sur la figure [2j le matériau est censé être sollicité en régime quasi- 
statique, et évoluer par une suite d'états d'équilibre. Or (comme il est indiqué au chapitre 9 du 
présentr trait 'e MIM), les simulations sont toujours rapides par rapport aux essais de laboratoire. 
Il est classique de caractériser l'évolution du système dans une telle simulation au moyen de 
diverses grandeurs liées au réseau des contacts telles que le nombre de coordination z (le nombre 
moyen de contacts par grain), la distribution des orientations des contacts, ou la distribution des 
forces. Toutes ces variables peuvent dépendre de l'écart à l'équilibre : un système mal équilibré 
comprendra en général moins de contacts, ceux qui portent de petites forces à l'équilibre pouvant 
facilement s'ouvrir sous l'effet d'une faible agitation résiduelle. Lorsque l'on cherche à savoir si le 
réseau des contacts est correctement déterminé, il est utile de tester la stabilité de l'équilibre au 
moyen du critère de positivité de la forme quadratique (l25j) . Cette opération se trouve facilitée 
si la contribution géométrique K*- 2 ) est négligeable (absence de mécanisme) et si le frottement 
n'est pas mobilisé. Or c'est effectivement ce qui arrive, comme on l'observe dans la pratique des 
simulations, du moins avec des disques ou des sphères : on constate que le nombre de contacts 
glissants tend à diminuer et à s'annuler à mesure que les écarts à l'équilibre se réduisent. Comme 
le système approche de son équilibre final, il est animé de différents modes de vibrations, dans 
lesquelles les forces de contact f oscillent de façon assez erratique, et, génériquement, l'arrêt de 
ces oscillations en un point de l'ensemble S des forces statiquement et plastiquement admissibles 
ne se produit pas sur sa frontière (là où l'inégalité de Coulomb est saturée dans au moins un 
contact). La positivité de la matrice de raideur, alors en bonne approximation symétrique, est 
assurée si les mouvements de mécanisme n'affectent pas la structure qui porte les forces. Pour 
choisir la tolérance sur l'équilibre on peut prendre une valeur de force telle que l'on obtienne avec 
des écarts plus faibles des matrices de raideur définies positives. On a ainsi observé, pour des 
assemblages de billes [10J qu'une tolérance de l'ordre de 10~ 4 -F/v, où Fn est la force de contact 
normale moyenne, est en général assez faible pour cela. Si on impose que le résultante des forces 
sur chaque grain soit inférieure à 10 _4 .F/v, que le moment résultant soit inférieur à 10~ F^d où 
d est le diamètre du grain, tandis que les contraintes extérieurement imposées sont équilibrées 
par les forces intérieures avec une erreur relative inférieure à 10 -4 , alors on vérifie que le réseau 
des contacts définit une matrice de raideur qui assure la stabilité su système. 

3.3 Hyperstaticité, hypostaticité, isostaticité et rareté des contacts 

Les réseaux de contact des assemblages granulaires à l'équilibre, dans le cas de grains 
circulaires (2D) ou sphériques, ne possèdent en général que des mouvements de mécanisme simples 
qu'il est facile d'éliminer lors de la construction des matrices de rigidité, quitte à diminuer la 
dimension Ni de l'espace des vitesses ou des petits déplacements. Ainsi, dans le cas d'une cellule 
périodique, les translations d'ensemble (mais pas les rotations) sont des mécanismes triviaux, 
au nombre de D. En général notons ko < D(D + l)/2 le nombre de mécanismes qui sont des 
mouvements d'ensemble de corps rigide. Une structure qui ne possède pas d'autre mouvement 
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de mécanisme est dite rigide (ou, plus précisément, rigide au premier ordres). 

D'autres mécanismes évidents sont les mouvements des grains flottants, c'est-à-dire qui 
ne transmettent aucune force. Dans les assemblages de disques ou de sphères faiblement polydis- 
persés, la proportion xq de grains flottants, selon la manière de préparer l'état d'équilibre peut 
varier de proche de zéro à plus de 15% |15[ IIP] , Avec une grande étendue granulo métrique |21| . 
le nombre de flottants parmi les grains de petite taille peut être beaucoup plus élevé. Tous les 
degrés de liberté des grains flottants sont des mécanismes, d'où k > uiXqN. 

Beaucoup d'études ont porté sur le cas des disques ou des sphères sans frottement [TJ [2j 
[22| Q3J dû]- Toutes les rotations sont alors des mécanismes, et on peut en fait les retirer de la 
liste des degrés de liberté car elles ne changent évidemment pas la géométrie du système et ne 
donnent lieu à aucun déplacement relatif normal dans les contacts. On peut retrancher 2DN C 
à la dimension des espaces de forces de contact ou de déplacements relatifs (en supprimant les 
composantes tangentielles), retrancher N à k, et remplacer Ni par = N — D(D — l)N/2 
(en supprimant les rotations), de sorte que (|10p devient 



Pour des objets non frottants de forme générale, les rotations ne sont plus a priori des mouvements 
de mécanisme, et on aura simplement 



tandis que pour des objets de révolution (3D) on a une rotation libre par grain et (|26p s'applique 
avec Jv/ 0) = Ni - N. 

La forme (j24]) de K*- 2 -* dans ce cas est une matrice symétrique négative, ce qui montre que 
tous les mécanismes dans lesquels il y a un déplacement relatif normal non nul dans au moins un 
contact conduisent à des instabilités. Si on a obtenu un état d'équilibre stable, alors on doit avoir 
k = ko + DNq dans (|27p . Nq étant le nombre de grains flottants. Dans un grand système, ceci 
entraîne une inégalité pour le nombre de coordination z* de l'ensemble des N* = N — Nq grains 
portant des forces à l'équilibre (z* = z/(l — xq)). En prenant Ne = 2z*N*, (127p . où k = ko est 
négligeable et h > donne z* > 6 pour des sphères (3D) et z* > 4 pour des disques (2D). En 
revanche, des grains de forme ellipsoïdale peuvent former sans frottement des assemblages stables 
avec des mouvements de mécanismes possibles non triviaux [23], la matrice K^ 2 ^ pouvant, selon 
la courbure des surfaces aux points de contact |11) . stabiliser ces mouvements ce qui autorise, 
selon (|10p . des valeurs z* < 12. 

Avec des grains frottants, l'unique mécanisme non trivial observé pour des billes sphé- 
riques est celui de la figure [3l qui met en mouvement les particules à deux contacts, le reste de 
l'assemblage restant fixe |10] , En chacun des deux contacts de la sphère mobile, il s'agit d'une 
combinaison de roulement et de pivotement, dans laquelle les forces de contact restent constantes, 
alors que les points de contact changent et décrivent sur la surface des grains une trajectoire cir- 
culaire. On vérifie que dans ce mouvement la contribution des contacts de la bille mobile au 
travail du second ordre est nul, le vecteur vitesse V est tel que K (1) • V = K^ 2 ) • V = 0. La 
stabilité ou l'instabilité ne se manifeste donc qu'au travers des effets d'une variation de moment 
subi de la part de la bille marquée 1 sur la figure par les billes marquées 2 et 3. Les simulations 
donnent une population faible (2 ou 3%) mais non négligeable de telles particules divalentes 
dans les assemblages faiblement coordonnés pjQ], l es forces qu'elles transmettent peuvent être 
importantes, et la proportion de mécanismes instables tend à disparaître quand la pression de 
confinement augmente. Si on trouve des grains de forme non sphérique avec deux contacts dans 

2. En général la propriété de rigidité au premier ordre est plus forte que la seule rigidité, qui est l'impos- 
sibilité de déformer la structure sans qu'il y ait déplacement relatif dans un contact [B]. 



N^' + h = N c + k. 



(26) 



N + h = N c + k 



(27) 
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FIGURE 3 - Mouvement de mécanisme d'une sphère divalente (numéro 1), la droite joignant ses 
deux points de contact avec les sphères 2 et 3, immobiles, est l'axe instantanje de rotation et 
porte les 2 forces de contact, (a) Equilibre des forces sur la sphère 1, dans le plan défini par 
les trois centres. L'équilibre est possible si a est inférieur à l'angle de frottement de contact, (b) 
Mouvement, vu de dessus (2 et 3 sont en position éclipsée). Le centre de la sphère mobile 1 décrit 
le cercle en pointillés fins autour de l'axe joignant les centres de 2 et 3. 

un assemblage 3D, on peut également leur associer un mouvement de mécanisme similaire, mais 
l'occurrence de telles configurations et leur stabilité ne semble pas encore avoir été répertoriées 
dans la littérature. 

En admettant l'absence de mécanismes non triviaux autres qu'associés aux grains di- 
valents, en proportion x\ parmi les N* = N(l — xq) grains qui transmettent des forces, la 
relation (|10p permet de minorer la coordinence z* = 2N C /N* du réseau des contacts actifs 



Pour le degré d'hyperstaticité, on dispose également de résultats pour les grains non frot- 
tants |24 [ l25 [ l26 | l27j. En effet, pour des configurations génériquement désordonnées (en pratique, 
pour tous les assemblages sauf les réseaux parfaitement ordonnés), le degré d'hyperstaticité des 
assemblages de grains non frottants est nul dans la limite rigide des faibles contraintes de confine- 
ment, ou des grandes raideurs de contact[f|. Une situation analogue familière est celle de la table 
à quatre pieds qui en général est bancale si ses contacts avec le sol sont rigides, parce que les 
incertitudes géométriques sur la forme des pieds comme sur les irrégularités du sol interdisent la 
configuration hyperstatique à quatre contacts. Cette propriété d'absence d'hyperstaticité est en 
fait de nature géométrique - le réseau des contacts ne peut pas supporter un système de forces 
normales auto-équilibrées - et sa validité est indépendante de la valeur effective du coefficient 
de frottement intergranulaire. Une conséquence immédiate en est une majoration du nombre de 
coordination z* de l'assemblage privé de ses grains flottants dans la limite rigide. A partir des 

3. C'est la limite oùk-} +oo, k étant le paramètre de raideur défini au chapitre 9. 




(28) 
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relations (|26p ou (|27p . on obtient en effet : 

z* < 2D (disques ou sphères, D = 2 ou 3) 

z* < D(D + l)/2 (grains quelconques, D = 2 ou 3) (29) 
z* < 10 (grains axisymétriques, D = 3) 

A la différence de (128p . ces inégalités restent vraies quel que soit le coefficient de frottement /â, 
mais seulement dans la limite des contacts rigides et indéformables. 

La structure des grains en contact portant des forces est isostatique quand elle est dé- 
pourvue d'hyperstaticité (h = 0) et d'hypostaticité, sauf les éventuels mouvements d'ensemble 
de corps rigide {k = ko). La matrice G, si on restreint l'espace des vitesses ou des petits déplace- 
ments en excluant de tels mouvements, est carrée et inversible. C'est effectivement le cas pour des 
billes ou des disques rigides non frottants (et non cohésifs), en présence de désordre générique. 
On a alors z* = 6 pour les billes, z* = 4 pour les disques, d'après f)28[) et (129p . Ces propriétés 
sont effectivement observées dans la limite rigide |28[ Q] |22l Q3J UU]. L'isostaticité est propre 
aux grains circulaires ou sphériques, des objets de forme différente pouvant s'assembler dans 
des configurations avec des mécanismes stables |23j . L'isostaticité a été exploitée pour mettre au 
point des méthodes de calcul quasi-statique d'assemblages de disques rigides se réarrangeant sous 
chargement variable sans aucun autre paramètre que géométrique |28j . L'absence d'hyperstati- 
cité s'applique plus généralement dans la limite rigide. Elle a des conséquences remarquables 
(les forces ne dépendent pas de la loi de contact) mais ne vaut pas pour les grains frottants, 
même dans la limite rigide |10| . sauf éventuellement pour certains procédés d'assemblage dans 
la limite \i — > oo (qui est une curiosité théorique). Il a parfois été suggéré qu'on pourrait avoir 
une sorte d'isostaticité « généralisée » pour des grains frottants à l'approche d'une rupture d'un 
assemblage, en donnant une coordonnée de moins au vecteur des forces de contact là où il y a 
glissement et donc égalité dans la condition de Coulomb. Cette prédiction est toutefois contredite 
par les observations numériques pQ. De plus, cette notion d'indétermination des forces prenant 
en compte le statut des contacts est d'un usage plus délicat, car elle dépend des forces appliquées 
et ne correspond à aucune propriété duale pour les déplacements. 

On notera que ces questions d'hyperstaticité ou d'hypostaticité ne prennent aucunement 
en compte les conditions inégalités qui portent sur les forces. Cependant, la majoration (I29p du 
nombre de contacts dans un assemblage granulaire désordonné entraîne une certaine limitation 
du degré d'hyperstaticité, qui est lié à certaines propriétés assez générales des matériaux granu- 
laires (comme la distribution des forces), et tend à restreindre l'influence de la loi de contact et 
à renforcer celle de la géométrie. Ainsi, en général, si l'on connaît avec précision les positions (et 
éventuellement les orientations) des grains, on peut calculer, avec les lois les plus habituelles, les 
forces de contact normales, qui sont élastiques et liées à la déflexion des contacts (qui appraît 
dans une simulation numérique comme une « interpénétration »). En revanche la force tangen- 
tielle dépend de l'histoire des sollicitations du contact et n'est pas directement déterminée par 
les positions actuelles. C'est pourquoi, dans la pratique des calculs par éleéments discrets, on ne 
sauvegarde pas seulement les positions mais aussi les forces à la fin du calcul. Cependant, on peut 
montrer que dans un état d'équilibre les forces tangentielles sont en général uniquement déter- 
minées pour les disques ou les sphères satisfaisant aux inégalités (|29p . Dans le cas de sphères, les 
équations d'équilibre sont au nombre de Ni = 6N* et contiennent z*N* coordonnées inconnues 
de forces tangentielles (c'est-à-dire 2 pour chacun des z*N*/2 contacts), et celles-ci sont donc, 
généralement, déterminées. Une conséquence en est que dans bien des évolutions quasi-statiques, 
l'oubli du dernier terme de (119p . qui viole le principe d'objectivité [TT] peut s'avérer inoffen- 
sif, la stabilité impliquant finalement un retour à l'équilibre avec les mêmes valeurs des forces 
tangentielles. 
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4 Applications à l'élasticité 



Nous rappelons ici quelques résultats relatifs aux propriétés élastiques évaluées numéri- 
quement pour des assemblages de billes pour lesquels l'élasticité des contacts obéit aux lois de 
Hertz-Mindlin-Deresiewicz [8j. Des comparaisons avec les résultats expérimentaux sont possibles, 
pour les valeurs des modules ainsi que pour l'extension du domaine de comportement (approxi- 
mativement) élastique. Elles donnent de bons résultats pourvu que les états internes du matériau 
dans l'expérience et dans la simulation soient proches, en particulier les nombres de coordination 
et les tenseurs de texture (orientation des contacts). Nous renvoyons à [29] 17] pour les références 
aux travaux expérimentaux sur les sables ou sur les billes et pour le détail des confrontations 
entre simulations et expériences. L'élasticité relie de petits incréments de contrainte à de très 
faibles déformations au voisinage immédiat d'un état de référence, et pourrait de ce fait être 
considérée comme peu pertinente pour le comportement mécanique des matériaux granulaires. 
Toutefois, les modules élastiques sont accessibles à l'expérience et fournissent indirectement des 
mesures non destructives de caractéristiques importantes des assemblages granulaires comme le 
nombre de coordination, qui, le plus souvent, échappent aux techniques d'observation directe [7J. 
La simulation de la réponse à de petites variations de chargement est en outre justifiée par l'étude 
de critères de localisation de la déformation [30], qui demandent une connaissance précise de la 
loi constitutive sous forme incrémentale. Dans cette partie notre propos est d'illustrer les apports 
de la méthode quasi-statique. Nous ne traitons pas des valeurs des modules, de leur dépendance 
par rapport aux contraintes et à la structure de l'assemblage |29[ 17], ni de leur prédiction, ou des 
propriétés quelque peu anormales des réseaux faiblement hyperstatiques |31[ HÏÏ] [7]. 

4.1 Calcul des modules élastiques 

De petits déplacements donneront des forces élastiques, comme on l'a vu au § 12.61 si on 
peut écrire une matrice de raideur symétrique définie positive. La cause majeure d'asymétrie 
dans la matrice de raideur est la mobilisation du frottement, Eq. (115p . Lorsque l'on a affaire à un 
assemblage bien équilibré on a en pratique ||Tjj|| < fiNij dans tous les contacts i-j (§ 13.2p . de 
sorte que l'on peut garder la forme symétrique et définie positive de chacun des blocs diagonaux de 
JC. L'élasticité suppose aussi que l'on puisse négliger la contribution géométrique non symétrique 
évaluée (pour des sphères) au § 12.51 ce qui est possible s'il n'y a pas de mécanisme. Lorsqu'il s'agit 
de calculer une réponse élastique, on ignore bien sûr les grains flottants. Restent les mécanismes 
associés aux grains divalents (figure [3]), qui appartiennent au noyau de et à celui de K^ 2 ), 
A moins d'exercer directement des efforts sur la particule divalente (marqué 1 sur la figure |3j), 
un incrément de chargement, en particulier une variation de contrainte globale, ne travaille pas 
dans ces mécanismes, que l'on peut éliminer soit en réduisant le nombre de degrés de liberté, 
soit en attribuant une raideur finie au mouvement libre, comme si, par exemple, le centre de la 
particule mobile était maintenue dans le plan de la figurera) par un ressort. Enfin, une dernière 
cause d'asymétrie de la matrice de raideur réside dans la sensibilité des champs de contraintes 
et de déformation dans la région du contact entre deux grains au trajet de chargement [9]. Il 
a été vérifié que ces subtilités de la loi de contact n'ont qu'un très faible impact sur le calcul 
des modules [7j. On peut donc à partir d'une configuration bien équilibrée construire la matrice 
de raideur élastique sous la forme K^ 1 ^ avec la forme élastique (|14p de la matrice des raideurs 
locales K. Après élimination éventuelle des mécanismes associés aux grains divalents, on résout 
le système linéaire 

g-U = AF ex \ (30) 

dont l'inconnue est le vecteur déplacement U, et l'incrément AF ext au second membre a des 
coordonnées nulles sauf celles qui expriment un incrément de contrainte (ou de force sur une 
paroi). Alternativement, on peut imposer des déformations de l'échantillon en manipulant les 
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conditions aux limites, ce qui revient à imposer certaines coordonnées de U. En partitionnant 
les degrés de liberté entre ceux qui sont imposés, coordonnées d'un vecteur u g de dimension 
n g <C N* et ceux qui sont laissés libres, coordonnées de U de dimension N? — n g , ainsi que la 
matrice K en blocs correspondants, (|30p prend la forme 
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(31) 



On trouve alors le système d'équations à résoudre pour U en prenant AF ext = dans (131 1) . où 
u g est connu : 

£-Û = -^-u 9 . (32) 

Que l'on utilise (|30p ou bien (|32p . on doit résoudre un système linéaire avec une matrice symé- 
trique positive définie et creuse. Pour ce faire, on dispose de diverses méthodes, dont celle du 
gradient conjugué |32j (éventuellement préconditionné), peu coûteuse en mémoire mais pouvant 
donner de longs calculs si la matrice de raideur est mal conditionnée, ou bien la factorisation de 
Cholesky, plus coûteuse pour un seul système linéaire, et pour laquelle il faut tenter de minimiser 
le stockage mémoire nécessaire en réordonnant au besoin les inconnues, mais qui finit par être 
avantageuse si on doit résoudre des systèmes avec la même matrice mais pour de nombreuses 
valeurs différentes du second membre. Les matrices de raideur des assemblages granulaires ayant 
une structure très similaire à celles que l'on rencontre dans les problèmes d'élasticité discréti- 
sés par éléments finis, on trouvera des indications utiles sur la résolution numérique dans la 
littérature plus vaste qui leur est consacrée. 

Bien entendu, il est possible de calculer des modules élastiques sans recourir à la matrice de 
raideur, en simulant par les méthodes dynamiques habituelles la réponse à de petits incréments 
de contrainte appliquée. Il faut ensuite s'assurer qu'il y a bien une réponse linéaire dans un 
certain intervalle de sollicitations, car des incréments de chargement trop faibles donnent des 
résultats affectés par les petits écarts à l'équilibre, et des incréments trop grands sortent du 
domaine élastique. Un tel calcul renseigne aussi sur le domaine élastique (ou approximativement 
élastique) . Le recours à la matrice de raideur |15[ [16] donne cependant un accès plus direct et 
plus rapide à l'ensemble des modules. 

4.2 Le domaine élastique 

La figure |4] montre les variations du déviateur des contraintes q = o\ — 03, normalisé par 
la contrainte latérale 03, et de la déformation volumique e„, près de l'état initial isotrope dans 
un essai de compression triaxiale de révolution simulé pour des billes de verre, en fonction de 
la déformation axiale e a = e\. L'indice 1 correspond ici à la direction principale majeure des 
contraintes. Les résultats sont montrés pour différentes valeurs de la contrainte isotrope initiale 
P = <7 2 = 03. On voit que l'élasticité au voisinage de l'état initial décrit en bonne approximation 
la relation entre incréments de déformation et de contrainte pour des déformations de l'ordre 
de 10 -5 , et que ce domaine approximativement élastique augmente avecP. Hors du domaine de 
validité de l'élasticité linéaire initiale, il est connu que la déformation est irréversible, comme le 
montre la figure [5] Au cours de la compression triaxiale, on voit que la pente de la courbe de 
déviateur décroît progressivement. Sous l'effet d'une décharge, on retrouve une pente plus élevée, 
proche du module élastique initial. En fait, la tangente à la courbe de décharge, on peut le vérifier, 
coïncide bien avec le module élastique correspondant, que l'on peut évaluer en faisant appel à 
la matrice de raideur pour l'état considéré. Ceci n'est possible que si la configuration est bien 
équilibrée. On doit donc, avant de calculer les modules par l'approche statique, laisser s'équilibrer 
l'état intermédiaire obtenu au cours de la simulation dynamique, en imposant des contraintes 
constantes plutôt qu'en poursuivant la déformation à vitesse contrôlée. Comme le montre la 
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FIGURE 4 - Déviateur des contraintes (à gauche), normalisé par la contrainte latérale 03 = P, et 
déformation volumique (à droite) fonctions de la déformation axiale dans la compression triaxiale 
simulée d'un assemblage de billes de verre pour 5 valeurs de P. Les points donnent les résultats 
du calcul par dynamique moléculaire, et les droites en lignes continues ont pour pentes E* et 
(1 — 2i/*), E* et v* étant le module d'Young et le coefficient de Poisson du matériau granulaire 
dans l'état initial isotrope. La flèche indique pour les courbes de e v le sens de P croissant. 




FIGURE 5 - Courbes q{e\)/a^ (axes de gauche) et e„(ei) (axes de droite), en compression triaxiale, 
simulée par une méthode de dynamique moléculaire classique, pour le système de la figure H] à 
P = 100 kPa (graphe de gauche), et pour un système similaire mais avec une coordinence initiale 
beaucoup plus grande (graphe de droite). Noter les échelles de déformation. Dans les deux cas, 
on a simulé les effets d'une décharge à partir de différents états atteints au cours de l'essai. Les 
lignes en pointillés fins représentent les résultats de la même simulation dans laquelle on a ignoré 
la création de nouveaux contacts. 

figure |6] cet équilibrage s'accompagne d'un léger fluage. Si on reprend par la suite la compression 
à taux de déformation axial imposé, on observe une remontée plus raide de la courbe de déviateur, 
dont la pente initiale coïncide avec le module élastique calculé par la matrice de raideur dans 
l'état d'équilibre atteint après fluage. Cette observation s'explique par la formation, lors de 
l'équilibrage, d'une structure plus stable avec les forces de contact strictement à l'intérieur du 
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FIGURE 6 - Équilibrage sous contraintes constantes, puis reprise de la compression triaxiale à 
k\ imposé à partir de différents états intermédiaires. Les symboles circulaires pleins montrent le 
début et la fin de chacun des intervalles de nuage qui correspondent à ces équilibrages (graphe 
de gauche). La reprise de l'essai à k\ imposé se caractérise par une réponse initiale très raide 
(courbes en trait gras) dont le graphe de droite, analogue de la figure 0] pour un état d'équilibre 
après fluage le long de la trajectoire de l'essai triaxial, montre qu'elle est élastique. 



cône de Coulomb, sous l'effet des vibrations du système autour de sa position finale, qui brouillent 
l'effet de polarisation des forces tangentielles du taux de déformation maintenue auparavant dans 
une direction constante. Une fois l'essai repris à k\ fixé, après une phase initial, les courbes de q et 
e v rejoignent celles de l'essai monotone et sans arrêt. Il est intéressant de noter que les modules 
élastiques se mesurent expérimentalement de manière similaire [33j : on applique de petits cycles 
de contraintes autour de la valeur à laquelle l'essai a été interrompu, ce qui provoque d'abord un 
certain fluage, puis les caractéristiques élastiques se déduisent de la forme finale stabilisée et peu 
dissipative du cycle de contraintes et déformations. En laboratoire c'est davantage la sollicitation 
appliquée qui est responsable du fluage que l'écart à l'équilibre@. 

Pour conclure, nous retiendrons de ces rappels rapides de résultats et d'observations numé- 
riques des propriétés élastiques des assemblages granulaires modèles que leur étude est facilitée et 
systématisée lorsque l'on a recours à l'approche quasi-statique fondée sur la matrice de raideur. 
Les conditions dans lesquelles on observe une réponse approximativement élastique et linéaire 
sont très similaires dans la simulation et dans les expériences : dans un intervalle du même ordre 
en très faible déformation suite au processus d'assemblage ; en décharge lors d'un essai triaxial : 
ou bien, sur de très courts intervalles, lors de la reprise de la déformation monotone dans la même 
direction suite à un arrêt d'un essai à é imposé et à un petit fluage provoqué par l'approche de 
l'équilibre ou par de faibles sollicitations cycliques. En dehors du régime approximativement 
élastique, on note que les déformations restent de l'ordre de la prédiction de l'élasticité linéaire 
d'autant plus longtemps que le réseau des contacts initiaux est bien connecté (figure [5]). 

4. On observe aussi un certain fluage dans les expériences mais seulement sur des durées beaucoup plus 
longues. On l'attribue au bruit ambiant. Voir le chapitre 9 pour une comparaison des valeurs de è entre les 
expériences et les simulations numériques. 
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5 Applications à la déformation anélastique 



5.1 Formulation du problème, propriétés générales 

Au-delà de la matrice de raideur élastique initiale, qui ne fournit que la tangente à l'ori- 
gine des courbes rhéologiques, voyons comment et dans quels cas l'approche quasi-statique peut 
prédire la déformation d'un échantillon granulaire soumis à un trajet de chargement donné. Pour 
fixer les idées nous traitons du cas simple de la compression biaxiale d'un échantillon bidimen- 
sionnel de disques, l'élasticité du contact faisant intervenir des raideurs Kn et Kt constantes. 
Nous présentons et commentons ici les résultats essentiels de nos propres travaux [U E| E] sur 
ce problème particulier, sachant que des éclairages très utiles sont fournis par les études de Mc- 
Namara et al. [5j [Ï3]. On traite le système dans l'HPP, on néglige K^ 2 - 1 , pour des contrainte 
appliquées de la forme a\ = P + q, oi = P, P étant la pression initiale isotrope. Le problème 
consiste donc à déterminer le vecteur déplacement U et les forces de contact f , tandis que le 
déviateur q augmente graduellement à partir de zéro. Le vecteur F ext ne contient que des va- 
leurs nulles, sauf les coordonnées qui correspondent aux contraintes, qu'elles s'expriment par des 
forces sur les parois (ai = F1/L2, 02 = F2/L1, voir la figure [2]) ou bien, dans le cas d'une cellule 
périodique, par (Aa a ) a= i t 2 qui doit satisfaire ©. En supposant que la trajectoire quasi-statique 
a été trouvée depuis l'état de départ jusquà une certaine valeur de q, pour laquelle le vecteur 
déplacement est U(çr), on cherche, pour un petit incrément de chargement AF ext correspondant 
à Aq le surcroît de déplacement AU, solution de : 

K(U(ç) , AU) • AU = AF ext . (33) 

Dans (|33p . la matrice de raideur K dépend de U et aussi de la direction de AU - on pourrait 

AU _ , . 

argumenter K par iT^jTj ■ On peut également écrire (j33|) en faisant apparaître les dérivées par 

rapport à q, qui joue le rôle d'un « temps cinématique », c'est-à-dire d'un paramètre le long de 
la trajectoire dans l'espace des configurations : 

KU,^)^-^, (M) 
— aq aq aq 

Comme on a une élasticité de contact unilatérale, la matrice de raideur est modifiée avant tout 
par la mobilisation du frottement, le bloc JÇ. prenant alors l'une ou l'autre des deux formes (fl~5j) . 
Elle est aussi affectée par l'ouverture des contacts - lorsque les grains i et j sont séparés, il faut 
bien sûr prendre jC,, = 0. Au cours de l'évolution quasi-statique des déplacements et des forces, 

un certain nombre N c de contacts deviennent « critiques » c'est-à-dire que le frottement y est 
complètement mobilisé. À une étape donnée du calcul, pour résoudre (|33p il faut déterminer le 
statut, glissant ou non glissant, de ces contacts critiques : a priori la matrice de raideur peut 
donc prendre 2 Nc formes différentes, et l'existence et l'unicité de la solution posent problème. 
Une propriété essentielle a été établie dans |13j : tant que l'ensemble des 2 Nc matrices de raideur 
satisfont le critère de stabilité, c'est-à-dire la positivité stricte de la forme quadratique (|25p . alors 
la solution existe et est unique. En pratique, plus la population de contacts glissants augmente, 
plus la stabilité est menacée. Lorsque certains statuts de contact pourraient donner A 2 W < 
avec certaines directions de <5U, on peut s'attendre à la manifestation d'une instabilité, un bruit 
arbitrairement faible pouvant solliciter le système dans la direction instable et déclencher une 
augmentation exponentielle de l'énergie cinétique. Il semble donc, grâce au résultat de |13| . que 
l'approche quasi-statique fournisse une solution unique aussi longtemps qu'elle est physiquement 
fondée. Cette conclusion peut être nuancée quelque peu cependant, car la démonstration de [13] 
ignore l'ouverture des contacts. À noter que ce phénomène n'est pas analogue au changement de 
statut : en effet, alors qu'un contact glissant peut devenir non-glissant à tout instant si la force 
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qu'il transmet pénètre à l'intérieur du cône de frottement, un contact ouvert ne se referme que 
si un incrément de déplacement fini vient combler l'interstice qui est apparu. L'ouverture des 
contacts est responsable (dans le cas linéaire) de la dépendance en U de K dans (|33p . alors que 
le changement de statut, de glissant à non glissant ou vice-versa, est la cause de sa dépendance 
dans la direction de l'incrément AU. 

Tant que |34]) peut déterminer une trajectoire comme une suite d'états d'équilibre, on 
a affaire à un réseau d'éléments rhéologiques (ressorts et patins frottants avec condition de 
Coulomb), la forme des grains n'intervient pas directement (pourvu que l'approximation HPP 
soit bonne, ce qui est effectivement le cas sauf pour des grains anormalement mous), sauf à 
conditionner initialement la géométrie du réseau de contacts. Il en résulte que les déformations, 
sous des contraintes données, seront pour une même géométrie inversement proportionnelles 
aux raideurs Kn, Kt- Nous qualifions ce comportement de régime I ou régime strictement 
quasi-statique. Il arrive aussi que l'évolution d'une collection de grains, qui reste proche de 
l'équilibre, se fasse par une succession de petits réarrangements |28[|2|[3]. qui sont déclenchés par 
des instabilités. Celles-ci ne donnent toutefois que des mouvements de faible amplitude, arrêtés 
par la fermeture de nouveaux contacts. Pour des systèmes de taille croissante, l'intervalle de 
contrainte (de q dans le cas de l'essai biaxial) est de plus en plus faible, ainsi que l'amplitude de la 
déformation qui accompagne le réarrangement du réseau j2], de sorte qu'à l'échelle macroscopique 
la déformation apparîtra comme graduelle et continue. Nous qualifions ce comportement, dans la 
limite des évolutions macroscopiques lentes, où l'inertie n'est plus pertinente^ de régime quasi- 
statique au sens large ou régime II - les déformations par réarrangement étant dites « de type 
II ». Nous donnons ci-dessous ( §5.3p quelques exemples d'utilisation de la méthode quasi-statique, 
ainsi que quelques illustrations des propriétés des régimes I et II. Auparavant, voyons comment 
on mène les calculs pour la résolution de ( 134]) . 

5.2 Un algorithme de calcul 

Pour résoudre numériquement (134p . on discrétise l'évolution du paramètre de chargement 
q en intervalle Aq et on utilise la forme incrémentale H33[) . Plaçons-nous pour une valeur donnée 
q, en supposant que le problème a été résolu pas à pas depuis q = 0, ce qui fournit les valeurs 
courantes de U et de f et cherchons à résoudre itérativement (|33p . La méthode de résolution que 
nous résumons ici se présente comme une recherche des incréments de forces de contact Af , tels 
que f + <5f soit plastiquement et statiquement admissibles (voir le § 13. ip . et de plus corresponde 
aux incréments de déplacement AU, qui doivent satisfaire la condition suivante. Définissant les 
déplacements relatifs élastiques par : 

AW E = (K?)- 1 ■ Af, 

avec la forme élastique de la matrice des raideurs de contact, et les déplacements relatifs plastiques 
AU?, par 

AU = T £ • AU = AÛ E + AÛ P , 

on doit avoir AZS = sauf pour les contacts glissants, où le déplacement relatif plastique est 
(avec nos conventions) positif^. Cette condition s'exprime aussi simplement en notant que l'on a 
f + !g-Û = V 
est défini sur 



f + /C E • U 



pour tout AU, où V désigne le projecteur sur le cône de Coulomb, qui 
a figure [7] En plasticité on dit qu'un tel projecteur définit la règle d'écoulement, 
c'est-à-dire le choix de la direction AU P . Un autre choix serait de projeter orthogonalement sur 
le cône, pour le produit scalaire qui corresponde à la norme ||f|| définie par 

iifii 2 = f • c^r 1 -f- (35) 



5. Voir au chapitre 9 comment le rôle de l'inertie est évalué par un nombre sans dimension. 

6. Le vecteur unitaire tangentiel est tel que Tij > 0. 
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FIGURE 7 - Projection V sur le cône de Coulomb suivant une règle d'écoulement non associée 
(graphe de gauche), ou associée (graphe de droite). 

correspondant au frottement de Coulomb usuel (graphe de gauche de la figure [7J est non associée. 
Une loi d'écoulement associée entraîne cette propriété remarquable que le calcul à la rupture 
fournit non seulement une condition nécessaire, mais aussi suffisante de stabilité (à condition 
toutefois que l'on puisse négliger l'influence de l'ouverture des contacts). Le réseau de contact 
continuera de supporter le chargement aussi longtemps qu'il existe des forces de contact qui 
soient à la fois statiquement et plastiquement admissibles. Pour le voir, il suffit de vérifier que la 
solution Af de (|33p minimise || Af|| 2 , avec la norme de (|35p . sous la contrainte f + Af € S. 

Pour résoudre itérativement f)33|) . on considère qu'on a affaire à un problème élastique 
corrigé par l'application de forces extérieures sur les grains concernés par les contacts glissants. 
Au début du calcul, lorsque l'indice d'itération j vaut zéro, on prend AU = K 1 AF ext , c'est- 
à-dire la solution élastique pour les incréments de déplacements. On évalue alors les forces de 
contact comme : 

fj = f + AUj, (36) 

ce qui donne des forces statiquement admissibles. Ces forces peuvent toutefois sortir du cône 
de Coulomb, aussi leur applique-t-on la projection V pour obtenir des forces plastiquement 
admissibles : 

ff A = ? H ■ (37) 
Le vecteur î? A n'est pas statiquement admissible, il n'équilibre pas F ext mais 

F ext + T^. (fPA-f.). 

On va par conséquent corriger le vecteur déplacement pour équilibrer ces forces (en admettant 
le problème élastique). Cette correction vaut 

Vj = g 1 • T £ • [f, - ff 4 ] = K -1 • [F ext - T £ • ff A ] , (38) 
et on incrémente AU : 

AU j+ i= AU j + Vj. (39) 

On peut alors remplacer j par j + 1, passant ainsi à l'itération suivante qui reprend au calcul 
des forces par (|36p . et on poursuit jusqu'à ce que Vj ou bien fj — fj soit négligeable. 

Cet algorithme, du point de vue des forces, revient à projeter alternativement sur le cône 
de Coulomb avec V, et sur l'espace affine des forces statiquement admissibles avec Q qui est une 
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projection orthogonale au sens de la norme (|35p . Les opérations H38[) et (|39p à l'étape j, suivies 
de (l36|) à l'étape j + 1 se traduisent par : 

f;+i = 2 RH ■ 

[Vérifions-le en définissant Q [f] pour un f quelconque. L'espace sur lequel il s'agit de 
projeter est l'ensemble des g tels que G • g = F ext . Il s'agit de décomposer f en Q [f] + f , 
avec f orthogonal au sens du produit scalaire associé à f)35|) au directeur de cet espace affine. 
En d'autres termes KT 1 ■ f doit être orthogonal au sens du produit scalaire ordinaire entre 
déplacements relatifs et forces de contact, à tout g tel que T G • g = 0, soit au noyau de T G. 
L'orthogonal du noyau de T G n'est autre que l'image de G, on a KT l • f = — G • V, où V est un 
certain vecteur de déplacements. Pour le déterminer, il suffit d'écrire que Q [f] est statiquement 
admissible. Appliquant l'opérateur G aux deux membres de l'égalité 

f = Q [f] - £ • <i • v, (40) 

on trouve 

v = K~ 1 (F ext - T G-f), 

et on en déduit Q [f] par (140p , Ce sont là exactement les opérations qui conduisent de à 

Le même algorithme, dans le cas associé, pour lequel V est un projecteur orthogonal au 
sens de la même norme (|35p que Q, permet de retrouver la propriété que si le chargement est 
supportable (c'est-à-dire si S 7^ 0), alors on trouvera effectivement une solution au problème 
élastoplastique, quelle que soit l'histoire du chargement. En effet il s'agit alors de projeter ortho- 
gonalement alternativement sur deux parties convexes fermées d'un espace de dimension finie, et 
la suite obtenue doit converger vers un élément de leur intersection S si elle n'est pas vide. 

Il n'en est pas de même avec la loi de glissement non-associée, car le déviateur maximal que 
l'on atteint est strictement inférieur au résultat du calcul associé, comme le montre la figure 0: le 
déviateur maximal atteint en régime I est voisin de 0, 8P dans le calcul non associé et supérieur 
à 1,3P dans le calcul associé. Dans la mise en œuvre de la méthode quasi-statique, on doit 
prendre en compte l'ouverture des contacts. En fait on peut le faire dans le cadre itératif de 
l'algorithme décrit ci-dessus, en gardant la matrice de raideur initiale et en faisant intervenir des 
forces auxiliaires pour corriger l'effet de l'annulation de certains blocs /Ç_. 

5.3 Illustration : régimes I et II en compression biaxiale 

La méthode quasi-statique permet de calculer l'évolution du système sous le trajet de char- 
gement biaxial jusqu"a une certaine valeur q\ du déviateur, qui borne l'intervalle de déformation 
de type I et vaut environ 0, 815P dans létude de [T], qui porte sur la compression triaxiale de 
systèmes de disques polydispersés assemblés initialement dans une configuration très dense (en 
fixant fi = pendant l'assemblage, cf. le chapitre 8 de cet ouvrage), avec une grande raideur de 
contact (K]\r/P = 10 5 ) et un coefficient de frottement \i = 0,25. q± ne diminue pas lorsque N 
augmente (à la différence de l'intervalle de stabilité des réseaux de contact entre grains rigides 
non frottants [28, 2j). Il semble, au contraire augmenter légèrement (les résultats sont compa- 
tibles avec un effet de taille finie en — (2, 12)/\/-W). Il est montré par ailleurs dans [TJ que q\ 
est indépendant de la raideur (si elle est assez grande, d'ordre 10 ) et du rapport Kt/Kj^. On 
observe par ailleurs qu'une proportion finie de contacts adopte le statut glissant (jusqu'à 20% 
dans [T]) et que 5 à 10% des contacts sont perdus. A noter que dans l'approche statique, aucune 
vibration « parasite » ne vient brouiller la distinction entre contacts glissants et non glissants. 

gi reste nettement inférieure au maximum de déviateur, ç/ max , d'où un régime de défor- 
mation par réarrangements (type II) pour q\ < q < ç max . La figure [9] illustre ces deux phases du 
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FIGURE 8 - Comparaison des régimes de chargement supporté en déformation strictement quasi- 
statique, pour la loi de contact habituelle, non associée (courbes marquées « n. a. »)et pour une 
loi de glissement associée (courbes marquées « a. »). 
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FIGURE 9 - Déviateur q en fonction de la déformation « axiale » dans un essai biaxial 2D 
simulé à contraintes contrôlées par petits paliers. En insert (noter la dilatation des échelles de 
déformation) on montre l'intervalle q < q\ en régime I, et on compare les calculs avec la méthode 
quasi-statique avec les résultats de simulation dynamique. 



comportement dans une compression biaxiale monotone. Elle représente le résultat de calculs par 
dynamique moléculaire à contrainte contrôlée, en imposant des pas de déviateur Aq = 10 -3 P, 
puis en attendant l'équilibre pour chaque valeur de q avant de l'incrémenter à nouveau. Lorsque 
le réseau initial des contacts reste stable, on a des déformations de type I d'ordre K^ 1 , comme 
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les déformations élastiques, et si faibles que la courbe se confond avec l'axe des ordonnées sur le 
graphe. En dilatant l'échelle des déformations (en insert sur la figure), on voit que dans ce régime 
les calculs en dynamique moléculaire et par la méthode quasi-statique sont en excellent accord. 
Au-delà de qi, on observe une courbe q{e\) en forme d'escalier. Dans les phases de stabilité 
(parties d'allure verticale) la déformation est à nouveau de type I et d'ordre KZ 1 . On a pu véri- 
fier qu'un calcul par la méthode quasi-statique était possible. Dans les phases de réarrangement 
(parties horizontales), le système se déforme par rupture du réseau des contacts jusqu'à ce qu'un 
nouveau réseau apparaisse et soit capable de supporter le déviateur appliqué. A la différence des 
déformations de type I, l'amplitude de ces événements de rupture n'est pas liée à la raideur des 
contacts. Les déformations (de type II) qui en résultent sont analogues aux déformations que l'on 
observe avec des modèles de grains rigides, comme en dynamique des contacts (voir le chapitre 
3 de ce traité). La sensibilité au niveau de raideur est d'ailleurs un moyen de détecter le type de 
déformation - voir à ce propos la discussion de l'influence du niveau de raideur sur le compor- 
tement quasi-statique au chapitre 9. Un autre moyen d'identifier, au moins approximativement, 
la nature (I ou II) des déformations est de tester jusqu'où il est possible, dans un calcul dyna- 
mique, de simuler, par exemple, un test biaxial ou triaxial, lorsque l'on ne crée aucun contact 
nouveau : on éprouve alors la stabilité du réseau initial. Le résultat de tels calculs est montré 
sur la figure [5] : l'intervalle de déviateur en régime I s'étend environ jusqu'à q = 0, 2o"3 pour le 
système initialement le moins coordonné, et jusqu'à q — 1, 103 pour l'assemblage de coordinence 
plus élevée. Il est naturel qu'un réseau de contacts mieux connecté soit capable de supporter un 
intervalle de contraintes macroscopiques plus étendu. 

À ce jour nous ne disposons pas d'analyse précise des mécanismes de rupture des réseaux 
de contact pour q = q\. C'est une perspective prometteuse dans l'étude fine des mécanismes de 
déformation des assemblages granulaires (à rapprocher d'autres matériaux amorphes). 

6 Conclusion 

Quoique loin de concurrencer les méthodes dynamiques, polyvalentes et d'emploi plus fa- 
cile, les approches quasi-statiques, fondées sur la construction de matrices de raideur, sont de 
précieux outils d'analyse des assemblages granulaires solides, à l'équilibre ou en déformation 
quasi-statique. D'un point de vue pratique, la construction de ces matrices fournit des moyens 
commodes pour juger de la stabilité des configurations d'équilibre et pour évaluer leurs pro- 
priétés élastiques. L'étude des matrices de rigidité et de raideur met en lumière les influences 
des différentes données géométriques et mécaniques et fournit d'utiles indications sur le nombre 
de coordination. Hors du petit domaine de réponse approximativement élastique, l'approche 
quasi-statique montre l'existence de deux régimes de comportement caractérisés par des origines 
physiques distinctes de la déformation macroscopique et des sensibilités différentes aux para- 
mètres micromécaniques. Encore assez embryonnaire, l'usage de l'approche quasi-statique et des 
matrices de raideur devrait trouver des applications fructueuses dans les études précises des mé- 
canismes de déformation des assemblages granulaires par instabilité, rupture et réarrangement à 
l'échelle microscopique. Comment se comporte la distinction entre régimes I et II dans la limite 
des grands systèmes et dans la limite des grains rigides? Avec des grains frottants de forme 
non sphérique, existe-t-il des mécanismes stables, sources de « modes mous » dans le spectre de 
vibration? Quelle est l'allure à grande échelle des champs de déplacements lors du déclenche- 
ment de la rupture? Comment le processus de déformation par rupture dépend-il de la forme 
des grains ? Telles sont certaines des questions assez fondamentales que le développement des 
méthodes quasi-statiques devrait permettre de clarifier. 
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